arXiv:math/0503219v3 [math.DG] 24 Feb 2006 


A discrete Laplace-Beltrami operator for simplicial 

surfaces 

Alexander 1. Bobenko Boris A. Springborn 
February 24, 2006 


Institut fill’ Mathematik, Technische Universitat Berlin, 

Strasse des 17. Juni 136, 10623 Berlin, Germany 

bobenkoSmath.tu-berlin.de 
springbSmath.tu-berlin.de 

Abstract 

We define a discrete Laplace-Beltrami operator for simplicial sur¬ 
faces (Definition 16). It depends only on the intrinsic geometry of the 
surface and its edge weights are positive. Our Laplace operator is sim¬ 
ilar to the well known finite-elements Laplacian (the so called “cotan 
formula”) except that it is based on the intrinsic Delaunay triangula¬ 
tion of the simplicial surface. This leads to new definitions of discrete 
harmonic functions, discrete mean curvature, and discrete minimal sur¬ 
faces. The definition of the discrete Laplace-Beltrami operator depends 
on the existence and uniqueness of Delaunay tessellations in piecewise 
flat surfaces. While the existence is known, we prove the uniqueness. 

Using Rippa’s Theorem we show that, as claimed, Musin’s harmonic 
index provides an optimality criterion for Delaunay triangulations, and 
this can be used to prove that the edge flipping algorithm terminates 
also in the setting of piecewise flat surfaces. 

Keywords: Laplace operator, Delaunay triangulation, Dirichlet energy, 
simplicial surfaces, discrete differential geometry 

1 Dirichlet energy of piecewise linear functions 

Let 5 be a simplicial surface in 3-dimensional Euclidean space, i.e. a geomet¬ 
ric simplicial complex in whose carrier S is a 2-dimensional submanifold. 
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in Berlin. 
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possibly with boundary. We assume S to be finite. Let V = {xi,... 

E, and F be the sets of vertices, edges and (triangular) faces of S. Let 
/ : S' ^ K, be a piecewise linear (PL) function on S (linear on each simplex 
of S). Then the gradient V/ is constant on each triangle. The Dirichlet 
energy E{f) = 5/5 ||V/|p is 

{xi,Xj)GE 

where the edge weights are 

_ J 5 (cot aij + cot aji) for interior edges 
I I boundary edges 

and aij, aji are the angle(s) opposite edge {xi,Xj) in the adjacent trian- 
gle(s) (see Figure 1). This formula was, it seems, hrst derived by Duf- 
fin [7], who considers triangulated planar regions. It follows (by sum¬ 
mation over the triangles) from the observation that the Dirichlet energy 
of a linear function on a triangle (xi,X 2 ,X 3 ) with angles ai, 02 , as is 
-^(fl(xi,X2,X3)) 4 XvigZ/SZ ~ /(^i+2)) • 



Figure 1: The a-angles of an internal edge. 

In analogy to the smooth case, the Laplace operator is defined as the 
gradient of the Dirichlet energy. (We identify the vector space of PL func¬ 
tions S' —> IR, with the the vector space IR,'^ of functions on the vertices.) By 
differentiating E(f) with respect to the value of / at a vertex Xi € V one 
obtains the “cotan formula” for the Laplace operator: 

Xj£V:{xi,Xj)£E 

Dziuk was the first to treat a finite element approach for the Laplace operator 
on simplicial surfaces, but without stating the cotan formula explicitely [9]. 
It seems to have been rediscovered by Pinkall and Polthier in their inves¬ 
tigation of discrete minimal surfaces [17], and turned out to be extremely 
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important in geometry processing where it found numerous applications, 
e.g. [6], [4] to name but two. In particular, harmonic parameterizations 
n : 14 ^ are used in computer graphics for texture mapping. The cotan- 
formula also forms the basis for a theory of discrete holomorphic functions 
and discrete Riemann surfaces [8] [15]. 

Two important disadvantages of this definition of a discrete Laplace 
operator are: 

1. The weights may be negative. The properties of the discrete Laplace 
operator with positive weights {wij > 0) are analogous to the properties 
of the classical Laplace-Beltrami operator on a surface with Riemannian 
metric. In particular the maximum principle holds. But some weights Wij 
may be negative, and this leads to unpleasant phenomena: The maximum 
principle does not hold. As a consequence, a vertex of a discrete minimal 
surface (as defined by Pinkall &: Polthier [17]) may not be contained in 
the convex hull of its neighbors. In texture mapping applications negative 
weights are undesirable because they lead to “flipped triangles”. In practice 
various tricks are used to avoid negative weights. 

2. The definition is not purely intrinsic. The classical Laplace-Beltrami 
operator is intrinsic to a Riemannian manifold: It depends only on the Rie¬ 
mannian metric. This is not the case with the discrete Laplace operator 
defined above. Two simplicial surfaces which are isometric but which are 
not triangulated in the same way give in general rise to different Laplace op¬ 
erators. As the simplest example, consider the two triangulations of a planar 
quadrilateral. They lead to different discrete Laplace operators. (Planarity 
is not what causes the problem since the quadrilateral may also be folded 
along either of its diagonals.) 

The key idea of this paper is that one can avoid both the above shortcom¬ 
ings by using the intrinsic Delaunay triangulation of the surface S to define 
the discrete Laplace operator (Definition 16) instead of the triangulation 
that comes from the simplicial complex S. 

2 Delaunay triangulations of piecewise flat sur¬ 
faces 

This section provides the necessary background on Delaunay tessellations 
of piecewise flat surfaces. We decided to give a detailed exposition because 
not all necessary proofs can be found elsewhere. 

The concept of a Delaunay triangulation in n-dimensional Euclidean 
space goes back to Delaunay [5]. Piecewise flat surfaces (Definition 1) were 
studied by (his student) Alexandrov [1] and more recently by Troyanov [20]. 
The idea of a Delaunay triangulation of a piecewise flat surface was ap¬ 
parently first considered by Rivin [19, Sec. 10]. The vertex set of the De¬ 
launay triangulation is assumed to contain the set of cone points of the 
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piecewise flat surface, so that the surface is flat away from the vertices. 
(This is very different from considering Delaunay triangulations in surfaces 
with Riemannian metric [13].) Rivin claims but does not prove an existence 
and uniqueness theorem for Delaunay triangulations in piecewise flat sur¬ 
faces. His proof that the edge flipping algorithm terminates is flawed (see 
the discussion after Proposition 12 below). A correct proof was given by 
Indermitte et al. [11]. (They seem to miss a small detail, a topological ob¬ 
struction to edge-flipability. See our proof of Proposition 11.) Furthermore, 
for our definition of the discrete Laplace-Beltrami operator we also need 
the uniqueness of the Delaunay tessellation, and this question has not been 
addressed properly. Rivin and Indermitte et al. define Delaunay triangula¬ 
tions by the local Delaunay criterion (see Definition 8), and infer existence 
via the edge flipping algorithm (see Proposition 12). To obtain uniqueness, 
we will define the Delaunay tessellation (whose faces are generically but not 
always triangular) via a global empty circle criterion (Definition 3). In a 
piecewise flat surface, the “empty circumcircles” are immersed empty disks 
(Definition 2) which may overlap themselves. Consequently, some work is 
required to show that this actually defines a (not necessarily regular) cell 
decomposition of the surface. Uniqueness, on the other hand, is immediate 
from this definition. We will also show that the local Delaunay criterion im¬ 
plies the global empty circumcircle condition. A Delaunay triangulation is 
obtained from the Delaunay tessellation by triangulating the non-triangular 
faces. It follows that a Delaunay triangulation, while in general not unique, 
differs from another Delaunay triangulation only by edges with vanishing 
cot-weights. This is important because it means that the discrete Laplace- 
Beltrami operator that will be defined in Section 3 depends only on the 
intrinsic geometry of the surface. 

Definition 1. A piecewise flat surface (PF surface) {S,d) is a 2-dimensional 
differential manifold S, possibly with boundary, equipped with a metric d 
which is flat except at isolated points, the cone points, where d has cone-like 
singularities. 

In other words, every interior point of a piecewise flat surface has a 
neighborhood which is isometric to either a neighborhood of the Euclidean 
plane or to a neighborhood of the apex of a Euclidean cone. The cone 
angle at the apex may be greater than 2tt. (Eor a more detailed definition 
of closed PE surfaces see Troyanov [20].) In this paper, we consider only 
compact surfaces and we require the boundary (if there is a boundary) to 
be piecewise geodesic. (The interior angle at a corner of the boundary may 
be greater than 27r.) 

A tessellation of a PE surface is a cell decomposition such that the faces 
are Euclidean polygons which are glued together along their edges. This 
implies that the cone points and the corners of the boundary are vertices of 
the cell decomposition. A triangulation is a tessellation where the faces are 
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triangles. For the following it is essential that we do not require tessella¬ 
tions (and in particular triangulations) to be a regular cell complexes, i.e. a 
glueing homomorphism may identify points on the boundary of a cell. For 
example, it is allowed that two edges of a face may be glued to each other; 
and an edge may connect a vertex with itself. A forteriori, we do not require 
a tessellation to be strongly regular, i.e. the intersection of two closed cells 
may not be a single closed cell. 

Remark. From the intrinsic point of view, the carrier S' of a simplicial surface 
S with the metric induced by the ambient Euclidean space is a piecewise 
flat surface. The simplicial surface S also provides S with a triangulation 
whose vertex set includes the cone points and the corners of the bound¬ 
ary. However, this triangulation is not intrinsically distinguished from other 
triangulations with the same vertex set. 

First, we will consider surfaces without boundary. To define the Delau¬ 
nay tessellation of a PF surface in terms of empty disks, we must allow an 
empty disk to overlap with itself: 

Definition 2 (immersed empty disk). Let [S, d) he a compact PF surface 
without boundary, and letV G S be a finite set of points which contains all 
cone points. An immersed empty disk is continuous map ip : D ^ S, where 
D is an open round disk in the Euclidean plane and D is its closure, such 
that the restriction ip\£) is an isometric immersion (i.e. every p £ D has a 
neighborhood which is mapped isometrically) and ip{D) n F = 0. 

Hence any points in ip~^iy) are contained in the boundary of D: (p~^iy) C 
dD. 

Definition 3 (Delaunay tessellation, no boundary). Let {S,d) be a 
compact PF surface without boundary, and let V C S be a finite set of 
points which contains all cone points. The Delaunay tessellation of {S, d) 
on the vertex set V is the cell decomposition with the following cells: A 
subset C C S is a closed cell of the Delaunay tessellation iff there exists an 
immersed empty disk ip ■. D ^ S such that (p~^{V) is non-empty and C is 
the image of the convex hull of ip~^{V): C = (^(conv (/?“^(F)). The cell- 
attaching map is ip\con-v ceZ/ is a 0-cell (vertex), 1-cell (edge), 
or 2-cell (face) if ip~^{V) contains one, two, or more points; respectively. 

Proposition 4. The Delaunay tessellation as defined above is really a tes¬ 
sellation of{S,d). 

Proof. Let us first remark that the vertex set of the Delaunay tessellation 
is obviously V. An edge e is a geodesic segment such that there exists an 
immersed empty disk tp : D ^ S with ip~^(V) containing exactly two points, 
ip~^iy) = {pi,P 2 \ such that e = ip{fii,P 2 \)-, where \pi,P 2 ] is the line segment 
joining pi and p 2 in D. That the vertices and edges form a 1-dimensional 
cell complex then follows from Lemma 5 below. 


5 


D 



Figure 2: The Intersecting Chord Theorem says that \\pi — q\\ ■ \\p 2 — g|| = 
||X 1 - Q'll • \\X2 - q\\. 

Then we have to show that the open faces are indeed homeomorphic 
to open disks. (A cell attaching map p : Iconv 93 - 1 ( 1 /) is a priori only an 
immersion of the interior of the domain.) This follows from Lemma 6 . 

It is comparatively easy to see that the cell attaching homeomorphism 
ip maps the boundary 9 conv{pi,p 2 ) • • • )Pn} into the 1 -skeleton; and that 
edges do not intersect open faces. We omit the details. 

Finally, Lemma 7 asserts that every point in S is contained in a closed 
cell. □ 

Lemma 5. The edges do not cross each other or themselves. 

Proof. Let e = p{\pi,P 2 ]) and e = <pi\pi,P 2 ]) be edges contained in the 
empty immersed disks p : D —> S and p : D ^ S, respectively, such 
that pi, p 2 are the only points in p~^{V) and pi, p 2 are the only points in 
p~^{y). Suppose e and e have an interior point in common: p{q) = p{q) 
with q G (pi,P 2 ) and q G (pi,P 2 )- We are going to show that 

p-^{e) = [pi,P2]- ( 1 ) 

Since p{pi),p{p 2 ) 0 p{D) the Intersecting Chord Theorem implies 

IIpi - 9II • \\P2 - q\\ < \\pi - q\\ ■ \\pi - q\\, 

where ||.|| denotes the Euclidean norm; see Figure 2. The opposite inequality 
follows equally. Hence p~^{p{pi)) and <^”^(<^(^ 2 )) must be contained in dD. 
Since pi and p 2 are the only points in p~^{y), this implies (1). □ 

Lemma 6. Letp ■. D ^ S he an immersed empty disk and suppose p~^{y) = 
{pi,P 2 , ■ ■ ■ )Pn} with n > 3. Let P = conv{pi,p 2 j • • • ,Pn}- Then the restric¬ 
tion y?|intp ofp to the interior of P is injective (and hence a homeomorphism 
int P —> (/?(int P)). 
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Figure 3; The dashed lines mark the boundaries of of the circular segments 
conv((dD) \ D) and conv((5i)) \ D). 


Proof. Suppose q € int P and q G D with q ^ q but ^p{q) = ^p{q). We will 
show that q 0 int P. 

Because ^p\d is an isometric immersion, there is a neighborhood q3 U C. 
D and an isometry T of the Euclidean plane such that T{q) = q, T{U) C D, 
and (p{T{x)) = (p{x) for all x G U. Let D = T{D) and cp = ipoT~^ : D —> S. 
Since (p and ip agree on the intersection D n D, there is a continuous map 
(p : Du D —f S such that ip\[) = p and (p\^ = (p. 

Now let Pi = T{pi) (i = 1,... n) and P = T{P) = conv{pi,p 2 , • • • ,Pn}- 
Then P and P have no common interior points: int P fl int P = 0. Indeed, 
since D and P are “empty”, {pi,p2, ■ ■ ■ ,Pn} C {dD)\D and {pi,P2, ■ ■ ■ ,Pn} C 
{dD) \ D. Hence P C conv((9P) \ D) and P C conv((9P) \ D). But the 
circular segments conv{{dD)\D) and conv((9P)\P) have no interior points 
in common. (See Figure 3.) 

Now g € int P implies q ^ int P. □ 

Lemma 7. Every point x G S is contained in a closed cell. 

Proof. Consider two immersed disks p:D^S,(p:D^S as equivalent if 
they differ only by a change of parameter, i.e. if = poT for some isometry 
of the Euclidean plane. The manifold of equivalence classes is parameterized 
by the set 

V = [{c,r) G S X Il>o I d{c, V) > r} 

of center/radius pairs. If we adjoin degenerated immersed empty disks with 
radius 0, we obtain a compact manifold with boundary, parameterized by 

V = {{c,r) G S X Il>o I d{c, F) > r}. 

For a point x G S the power function with respect to x is the continuous 
function 

pow,j, : P ^ H, pow,j,(c, r) = (d(c, x)) - r^. 
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Figure 4: If p is not contained in the convex hull of p~^(y) = {pi,... ,Pn}j 
then there is another immersed empty disk with smaller power. (See proof 
of Lemma 7.) 


If <y 9 : Z) —> S' is an immersed empty disk with center c and radius r, 
then P 0 W 3 , (c,r) is smaller than, equal to, or greater than zero depending 
on whether x € x G dip{D), or x € S \ p{D)- 

Let X € S. We have to show that x is contained in some closed cell. 
If X G IL this is clear, because the points in V are 0-cells. So assume 
x G S \ V. Since pow^, is continuous on the compact set V, there is a 
(cminj ^’min) £ ^ where pow^, attains its minimum. Since x ^V, there is an 
empty disk containing x and hence poWj-fcmin. < 0. Let p ■. D ^ S he 
an immersed empty disk with center Cmin and radius rmin; i.e. ZZ is a disk in 
with center m G IR,^ and radius rmin and <p(m) = Cmin- There is a p G ZZ 
with (p(p) = X and ||p — m\\ = d{x, C min ). We show by contradiction that 
p G conv(p“^(I/). 

Suppose the opposite is true; p 0 conv(p“^(F). Then there exists a 
closed half-space ZZ C with C ZZ but p ^ H. In that case 

there exists another immersed empty disk (p : D ^ M with D\D C ini H, 
D\D and (see Figure 4). Let fh and f be center 

and radius of D. Then 

Up — m|p — f^ < Up — m||^ — r^. 

(To see this, note that q 1 —> (||g — m|p — r^) — (||g — m|p — r^) is an affine 
linear function —> E, which vanishes on dH and is positive on int ZZ, 

and negative on E^ \ ZZ.) This implies poWj,(c, r) < pow,^(c m;n . Xmin), where 
c = (p{rh). This contradicts the assumption that pow^, attains its minimum 
on (Cminj ^min)' D 

Delaunay tessellations of PF surfaces with boundary. Now let (S', d) 
be a compact PF manifold with piecewise geodesic boundary. Let V C S 
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be a finite set of points which contains all cone points and all corners of the 
boundary. The boundary is then the union of geodesic segments connecting 
points in V but containing no points of V in their interior. Each connected 
component of the boundary is a closed geodesic polygon with vertices in V. 
To each such boundary polygon glue a PF surface obtained by cyclicly gluing 
together the appropriate number of isosceles triangles with appropriate base 
lengths and legs of length R > 0. Each of these caps contains a special point 
where the triangle apices are glued together. It is in general a cone point. 
The result of closing all wholes with such caps is a closed PF surface (S, d). 
Let V be the union of V and the set of special points of the caps. If R is 
chosen large enough, then the isosceles triangles in the caps will be faces of 
the Delaunay tessellation of (S, d) with respect to V. (This is so because if 
R is large enough, the immersed circumdisks intersect S in lunes which are 
so small that they are empty.) The Delaunay tessellation of the bounded 
surface (5, d) with respect to V is defined to be the cell complex obtained 
by removing these triangles. 

Delaunay triangulations. A Delaunay triangulation is a triangulation 
obtained from a Delaunay tessellation by triangulating all non-triangular 
faces. A Delaunay triangulation is characterized by the empty circumcircle 
property: A tessellation of {S, d) on the vertex set E is a Delaunay triangula¬ 
tion iff for each face / there exists an immersed empty disk ip : D ^ S such 
that there are three points pi,P 2 ,P 3 S with / = p{co\\\{pi,p 2 ,pz})■ 

(But there may be more than three points in ip~^iy).) 

In Proposition 10 we give a more local characterization of Delaunay 
triangulations. 

Definition 8. Let T be a geodesic triangulation of {S, d) with vertex set 
V, and let e be an interior edge ofT. Since all faces ofT are isometric to 
Euclidean triangles we can isometrically unfold to the plane the two triangles 
ofT that are adjacent to e. We say that e is locally Delaunay if the vertices 
of the two unfolded triangles are not contained inside the circumcircles of 
these triangles. 

For our investigation of discrete Laplace operators we will need the fol¬ 
lowing characterization of Delaunay edges. 

Lemma 9. An interior edge e of a triangulation T of a piecewise flat surface 
(5, d) is locally Delaunay if and only if the sum of the angles opposite e in 
the adjacent triangles does not exceed vr. 

This follows immediately from the fact that opposite angles in a circular 
quadrilateral sum to vr. 

Clearly all interior edges of a Delaunay triangulation are locally Delau¬ 
nay. This property actually characterizes Delaunay triangulations: 
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Figure 5: Left: The layed out triangles. Right: The power line of Cj and q+i 
is Since the edge is assumed to be Delaunay, Paix) > Pa+iix) 

in 


Proposition 10. A triangulation T of a piecewise flat surface {S, d) is 
a Delaunay triangulation if and only if all interior edges of T are locally 
Delaunay. 

The following proof is an adaptation of Delaunay’s original argument [5] 
for Delaunay triangulations in IR,"'. (See also Edelsbrunner [10, p. 8] for 
a more easily available modern exposition.) Alternatively, one could also 
adapt the argument of Aurenhammer and Klein [2] for Delaunay triangula¬ 
tions in the plane. 

Proof. If S' is a manifold with boundary, construct a closed PF surface by 
glueing piecewise flat disks to the boundary components in the manner de¬ 
scribed above. If R is chosen large enough, all edges will be locally Delaunay. 
It remains to prove the Proposition for closed PF surfaces. 

Suppose that all interior edges of the triangulation T of the closed PF 
surface (S, d) are locally Delaunay. We want to show that T is a Delaunay 
triangulation. To this end we will show that the empty circumcircle property 
holds. Let to be a triangle of T. Starting with to we develop a part of T 
in the Euclidean plane (see Eigure 5, left). Begin with a triangle tq in 
the Euclidean plane that is congruent to to- Let cq be the circumcircle of 
tq. Next, lay out congruent copies of the triangles neighboring to. This 
introduces new vertices in the plane, which are on or outside co because 
the edges of to are Delaunay by assumption. Keep laying out neighboring 
triangles in the plane at free edges but only if the free edges intersect co. 
(Different layed out triangles may correspond to the same triangle in (S, d).) 
Each new triangle introduces a new vertex and we will show that they do 
not lie inside co. Hence, when the layout process stops (when all free edges 
do not intersect co), the triangles simply cover the inside of cq. It follows 
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that to has an empty circumdrcle in (5, d). 

Let To,Ti,... ,Tn, be a sequence of layed out triangles such that Tj and 
Tj+i share an edge Let Xj+i be the vertex opposite in Tj+i. 

Assuming that Xi is not inside cq for all z < n we will show that the same 
holds for Xn- Let be the straight line containing and let 

be the half space bounded by and containing Tj+i. Then 

n Do) D {Hi ^2 n Do) D ... D {Hn-i,n n Do), 

where Do is the open disk bounded by cq. Hence it remains to consider the 
case where Xn € for all i = 0 ,..., n — 1, because otherwise Xn ^ Dq. 

Now consider the power of a point x G with respect to a circle c C 
with center Xc and radius r as a function of x: 

Pc{x) = \\x — XclP — 

It is positive, zero, or negative if x lies outside, on, or inside c, respectively. 
The power line of two different circles c and c' is the locus of points x with 
Pc{x) = Pc'{x)- It is a straight line because pdx) — Pc'd) is linear in x. 
The power line of two intersecting circles is the line through the intersection 
points. Let q be the circumcircle of Tj. Either Ci = Cj+i or the power line 
of Ci and c^+i is gi^i+i and 

Hi^i+i = {x :pcdx) >Pci+dx)}. (2) 

(See Figure 5, right.) Indeed, pci+i(3;i+i) = 0 and since the edge is 

locally Delaunay by assumption, pc^Xi+i) > 0. Hence 

PcoiXn) > PciiXn) > ■ ■ ■> Pcr,{Xn) = 0, 

and therefore Xn lies on or outside cq. This concludes the proof. □ 

The edge flipping algorithm may be used to construct a Delaunay trian¬ 
gulation of a piecewise flat snrface {S, d) with marked points H C S': 

1. Start with any triangulation T of (S, d) with vertex set V. 

2. If all interior edges of T are locally Delaunay, stop. 

3. Otherwise there is an interior edge e of T which is not locally Delaunay. 
Perform an intrinsic edge flip'. Replace e by the other diagonal of the 
quadrilateral formed by the two triangles adjacent to e. Go to Step 2. 

The following two propositions show that this is indeed an algorithm. 

Proposition 11. If an edge is not locally Delaunay, then it can be flipped. 
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Proof. In an edge is flipable iff the two adjacent triangles form a convex 
quadrilateral. In a PF surface, there is an additional topological obstruction 
to flipability: Since the triangulation may not be regular, an edge may be 
adjacent to the same triangle on both sides. So suppose the edge e is not 
locally Delaunay, i.e. the sum of opposite angles exceeds vr. Then there are 
two different triangles adjacent to e, because the sum of all angles in one 
triangle is vr. These two triangles form a Euclidean quadrilateral (possibly 
with some of the boundary edges identified with each other), which is convex 
by the usual argument. Hence e can be flipped. □ 

Proposition 12 (Indermitte et al. [11]). The edge flipping algorithm 
terminates after a finite number of steps. 

Together with Proposition 10 this implies that the edge hipping algo¬ 
rithm produces a Delaunay triangulation (in the global empty-circumcircle 
sense). To prove Proposition 12, one has to show that the algorithm cannot 
loop infinitely. In the setting of planar Delaunay triangulations, it is enough 
to define a suitable real valued function on the set of triangulations on the 
given vertices which decreases (or increases) with each edge flip. Because 
this set of triangulations is finite, the algorithm has to terminate. As a 
further consequence such a function attains its minimal (or maximal) value 
on the Delaunay triangulations. Several such functions are known, see for 
example Lambert [12], Musin [16], Rivin [19, Sec. 10], and the survey arti¬ 
cle [2]. When we consider Delaunay triangulations of PF surfaces, however, 
the set of triangulations on the marked points may be infinite. (The fact 
that the number of combinatorial types of triangulations is finite [19] is 
not sufficient to make the argument.) For example, the surface of a cube 
has infinitely many geodesic triangulations on the eight vertices. To prove 
Proposition 12 by means of a function which decreases with every flip, one 
has to show that it has the following additional property. 

Definition 13. Let T be the set of triangulations of a PF surface on a given 
set of marked points and let f : T — R. We say that f is proper if for any 
Af G IR, the number of triangulations T with f{T) < M is finite. 

In their proof of Proposition 12, Indermitte et al. [11] use the sum of 
squared circumcircle radii as proper function which decreases with every flip. 
(For a proof that it decreases with every flip, they refer to an unpublished 
PhD thesis. However, see Musin [16], his Theorem 3 and the following 
Lemma, for a hint on how to prove this.) Another possible choice is Musin’s 
harmonic index. Below we show that it is proper and decreases with every 
flip. The latter fact we deduce from Rippa’s Theorem [18], which is also 
of independent interest in connection with the Laplace-Beltrami operator. 
Rippa’s proof holds without change also for piecewise flat surfaces. 
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Rippa’s Theorem. Let (5, d) be a piecewise flat surface and let V C S 
be a set of marked points which contains the cone points and the corners 
of the boundary. Let f : V ^ M, be a function on the marked points. For 
each triangulation T of {S, d) with vertex set V let fx ■ S ^ M, be the PL 
interpolation of f that is linear on the faces of the triangulation T. 

Suppose Ti is a triangulation with an interior edge e and T 2 is obtained 
from Ti by flipping e. If the edge e is a Delaunay edge after the flip, i.e. in 
T 2 , then 

E{fT,)>E{fTfl, 

where E denotes the Dirichlet energy as in Section 1. Equality holds only if 
/ti = /t 2 or if e was also a Delaunay edge in Ti. 

As a consequence, the minimum of the Dirichlet energy among all pos¬ 
sible triangulations is attained on the Delaunay triangulations {S,d): 

min [ I V/t P= [ \ V/r^ p, 

^ Js Js 

where Td is any Delaunay triangulation. 

Moreover, for generic f : V ^ M., this property of Delaunay triangula¬ 
tions is characteristic. 


Rippa’s proof is based on the following comparison formula for the 
Dirichlet energies of two possible triangulations of a quadrilateral 


Eifrfl-EifTfl 


(/i - /2)^ in + r3){r2 -I rj) 
2 sin 0 


(nra - r2r4). 


Here Ti and T 2 are the two triangulations of the convex quadrilateral Q = 
(xi,X 2 ,^ 3 , 2 : 4 ) obtained by addition of the diagonals (xijXs) and (x 2 ,X 4 ) 
respectively, /i and /2 are the values of fxi and /tj at the intersection 
point xq of the diagonals, ri,..., r 4 are the distances from xq to the vertices 
xi,..., X 4 of the quadrilateral and 9 is the intersection angle of the diagonals. 

For a Euclidean triangle t with sides a, b, c and area A, Musin [16] defines 
the harmonic index as 


hrm{t) 


A 


The harmonic index of a triangulation T with face set F is the sum of the 
harmonic indices of all triangles: 


hrm{T) = hrm{t). 
teF 


Proposition 14. The harmonic index is proper. 
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Proof. In a PF surface, there may be infinitely many geodesic lines connect¬ 
ing two points, but for any L G IR, only a finite number of them have length 
< L [11]. Hence there are only a finite number of triangulations of a PF 
surface with marked points such that all edges are not longer than L. Now 
the Proposition follows from the inequality 

hrm{T) > 

^tot 

where is the largest length of an edge of T and Atot is the total 

area of the PF surface. Indeed, if hrm{T) < M, then lma.x{T) < MAtot, 
and there are only finitely many triangulations satisfying this bound on edge 
lengths. □ 

The following theorem was stated by Musin without proof [16]. 

Theorem 15. With the notation and under the conditions of Rippa’s The¬ 
orem 

hrm{Ti) > hrm{T 2 ) 

and equality holds only if e is a Delaunay edge in Ti as well. This implies 
that the harmonic index is minimal for a Delaunay triangulation To (and 
hence for all of them): 


imn hrm{T) = hrm{To). 

Proof. The harmonic index of a triangle is 

hrm{t) = 4(cot a + cot (3 + cot 7 ), 

where a, (3, 7 are the angles of the triangle. (Because the area is H = 
^aha = \hhh = \chc, where ha,hij,hc are the heights of the triangle; and 
a/ha = cot (3 + cot 7 , etc.) 

For a triangulation T of {S, d) with vertex set V and x G P let dx,T ■ 
5 —> IR, be the function that is linear on the triangles of T, equal to 1 at x, 
and equal to 0 at all other marked points in V. Then 

^ F;(4,r) = ^ ^ cot a = i hrm{T). 

x£V angles q in T 

Hence the theorem follows from Rippa’s Theorem. □ 

3 The discrete Laplace-Beltrami operator and dis¬ 
crete harmonic functions 

We are now in a position to define the discrete Laplace operator on a sim- 
plicial surface in an intrinsic way. 
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Definition 16. Let S be a simplicial surface with vertex set V and let S be 
its carrier, which is a piecewise flat surface. The discrete Laplace-Beltrami 
operator A of a simplicial surface S is defined as follows. For a function 
f : V ^ on the vertices, the value of A/ : V ^ at Xi £ V is 

^f{xi)= ^ixi,Xj){f{xi) - f{xj)), (3) 

Xj£V-.{xi,Xj)&Eo 

where Ed is the edge set of a Delaunay triangulation of S and the weights 
are given by 

, . f i (cot 041+ cot a,-j) for interior edges 

= S 1 , r 1 1 7 ■ (4) 

I 2 boundary edges 

Here aij (and aji for interior edges) are the angles opposite the edge (xi,Xj) 
in the adjacent triangles of the Delaunay triangulation (see Figure 1). 

The discrete Dirichlet energy of / is 

= X] '^{Xi,Xj){f{Xi) - f{Xj)f. 

{xi,Xj)^Eo 


Due to Lemma 9 and the formula cot a. + cot B = the dis- 

Crete Laplace operator has non-negative weights. The edges with vanishing 
weights are diagonals of non-triangular cells of the Delaunay tessellation. 
Erasing such edges in (3) we obtain a discrete Laplace operator on the De¬ 
launay tessellation of S with positive weights n{x,Xi). Moreover, this prop¬ 
erty is characteristic for Delaunay triangulations: Consider a piecewise flat 
surface {S, d) with a triangulation T. Denote by Ay the Laplace operator 
of the triangulation T: it is given by the same formula (3) with the weights 
ut determined by the triangulation T by the same formulas (4) as for the 
Delaunay triangulation. The following observation is elementary. 

Proposition 17. The Laplace operator At of the triangulation T has non¬ 
negative weights ut if and only if the triangulation T is Delaunay. 

Laplace operators with positive weights on graphs possess properties 
analogous to the smooth theory. 

Definition 18. A discrete function f : V ^ on a simplicial surface is 
harmonic if Af{x) = 0 for all interior vertices x. 

Discrete harmonic functions satisfy the maximum principle: A real val¬ 
ued harmonic function attains its maximum on the boundary. This implies: 

Proposition 19. For each interior vertex x, the value f{x) of a harmonic 
function f : V —>■ lies in the convex hull of the values f{xi) on its 
neighbors. 
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We conclude this section with some standard facts regarding boundary 
value problems for the discrete Laplace operator. Let us fix a subset Vq <ZV 
of vertices-which may but need not be the set of boundary vertices-and a 
function g : Vg ^ IR,. The problem of finding a function that satisfies 
A/(x) = 0 for X G y \ Va and f{x) = g{x) for x G Va is called a Dirichlet 
boundary value problem. The problem of finding a function that satisfies 
A/(x) = 0 for X G y \ Va and A/(x) = g{x) for x G Va is a Neumann 
boundary value problem. 

Theorem 20. For arbitrary Vg and g{x) the solutions of the Dirichlet and 
Neumann boundary value problems exist and are unique (up to an additive 
eonstant ifVg = $ and in the Neumann case). The solution minimizes the 
Dirichlet energy £d. 

Indeed, solutions of the Dirichlet and Neumann boundary value problems 
are critical points of the Dirichlet energy. Since the energy is a positive 
definite quadratic form, its only critical point is the global minimum. 

More complicated boundary conditions, such as so called “natural bound¬ 
ary conditions” (see Desbrun et al. [6]), are also intensively used in geome¬ 
try processing. The corresponding existence and uniqueness results are still 
missing. 

4 Discrete mean curvature and minimal surfaces 

In this section we adapt the definitions of the mean curvature vector for 
simplicial surfaces and minimal surfaces originally suggested by Pinkall and 
Polthier [17] to the discrete Laplace operator that we introduced in Section 3. 

For a smooth immersed surface / : IR,^ D 1/ —> IR^ the mean curvature 
vector is given by the formula H = A/, where A is the Laplace-Beltrami 
operator of the surface. For a simplicial surface we define the mean curvature 
vector by the same formula, following [17], but we use a different Laplace 
operator: 

Definition 21. Let S be a simplieial surface with carrier S. The discrete 
mean curvature vector at a vertex x is 

nix) = A/(x), 

where f : S ^ 'E? is the restriction of the identity map on IR^ to S, and A 
is the discrete Laplace-Beltrami operator of Definition 16. 

The discrete mean curvature vector 7f(x) of a simplicial surface corre¬ 
sponds to the integral of the mean curvature vector of a smooth surface over 
a neighborhood of the point x. (Note when we scale the simplicial surface, 
n varies as the integral of the mean curvature over some domain.) This 
suggests the following alternative definition: 
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Definition 22. In the setup of Definition 21, let C{x) be a Voronoi cell 
of the vertex x of a simplicial surface. The discrete mean curvature vector 
density at x is defined by 


H{x) 


Af{x) 

A{Cix)y 


where A is the discrete Laplace-Beltrami operator of Definition 16 and A[C{x)) 
is the area of the Voronoi cell C{x). 


This definition is similar to the definition of mean curvature suggested 
by Meyer, Desbrun, Schroder and Barr [14]. Again, the difference is that we 
contend that one should use the discrete Laplace-Beltrami operator. 

Definition 23 (Wide definition of simplicial minimal surfaces). A 

simplicial surface is called minimal if its mean curvature vector vanishes 
identically. 

Since the embedding / : 5 ^ of a simplicial minimal surface is 
harmonic, Proposition 19 implies 


Proposition 24. Every interior vertex of a simplicial minimal surface lies 
in the convex hull of its neighbors. 


The following stricter definition is also natural. 

Definition 25 (Narrow definition of simplicial minimal surfaces). 

A simplicial surface S is called minimal if its mean curvature vector van¬ 
ishes identically and the intrinsic Delaunay triangulation of the carrier of S 
coincides with the triangulation induced by the simplicial complex S. 


Such a simplicial minimal surface is a critical point of the area functional 
under variations of the vertex positions [17]. We would like to note that 
there exists also a non-linear theory of discrete minimal surfaces based on 
the theory of circle patterns [3]. 

The mean curvature flow for simplicial surfaces is given by the equation 

|(x)=//(x). 

This flow may change the Delaunay triangulation of the surface. At some 
moment two Delaunay circles coincide and the diagonal of the corresponding 
quadrilateral flips. However, since the weights of the diagonals vanish at the 
flip moment, the discrete Laplace-Beltrami operator and mean curvature 
vectors are continuous functions of time t. 

For the numerical computation of discrete minimal surfaces one should 
use the algorithm of [17] with an extra step to adapt the Delaunay triangu¬ 
lation: Start with a simplicial surface Sq which respects the given boundary 
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conditions. Calculate the Delaunay triangulation of the carrier Sq and the 
weights u. Find an / : 5o —> IR,^ which respects the boundary conditions 
and minimizes the Dirichlet energy (see Dehnition 16). You may start with 
the embedding of Sq as initial guess. You obtain a new simplicial surface 5i 
which is combinatorially equivalent to Sq but geometrically different. Calcu¬ 
late the Delaunay triangulation and weights of Si and hnd an / : 5i —> IR,^ 
that minimizes the Dirichlet energy. Iterate. 
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